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Abstract 

The Lagrangian derivatives of finite-time Lyapunov exponents and the corre- 
sponding characteristic directions are shown to satisfy time-asymptotic differential 
constraints in chaotic flows. The constraints are valid for any metric tensor, and 
are realised with exponential accuracy in time. Some of these constraints were de- 
rived previously for chaotic systems on low-dimensional Euclidean spaces, by requir- 
ing that the Riemann curvature tensor vanish in Lagrangian coordinates. The new 
derivation applies in any number of dimensions, predicts the number of constraints 
for a given flow, and provides a rigorous convergence rate of the constraints. 



1 Introduction 



A deep understanding of chaotic systems originates from the study of time- 
asymptotic dynamics [1]. The Lyapunov exponents (or characteristic expo- 
nents) are the most well-known instance; they describe the average rate of ex- 
ponential separation of neighbouring trajectories, in the infinite-time limit [2- 
4]. The Lyapunov exponents are often used as a criterion for chaos: on a 
bounded domain, at least one positive exponent is required. They have proven 
quite rich, for example in characterising the fractal dimension of attractors in 
dissipative systems [5,6]. 



The convergence of the Lyapunov exponents is extremely slow, typically log- 
arithmic in time. This is due to the fact that they are inherently nonlocal 
and require the fiducial trajectory (the trajectory along which the exponents 
are evaluated) to explore the attractor (or invariant region, for nonattracting 
systems) very thoroughly, as the greatest stretching often comes from rare 
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excursions near hyperbolic fixed points. It is thus debatable whether the ex- 
ponents are relevant to the finite-time dynamics of a single trajectory, or to 
an ensemble of tightly bunched trajectories, since a laboratory or numerical 
experiment will usually have ended long before the exponents have converged. 

Hence the shift of emphasis to the finite-time Lyapunov exponents (FTLEs). 
These are the same as the Lyapunov exponents (i.e., the infinite-time Lya- 
punov exponents) but averaged over a finite time, as the name implies. On 
short to moderate time scales, the FTLEs do not show convergence at all 
(that is, the average rate of separation is not exponential), but they are still 
connected to a rapid, quasi-exponential average separation of neighbouring 
trajectories. 

The rapid separation of trajectories leads to a quasi-exponential convergence 
of the characteristic directions, the directions in the space of initial conditions 
(Lagrangian coordinates, which label fluid elements) exhibiting a given char- 
acteristic rate of separation. Thus, even though the FTLEs have not usually 
converged, on moderate time scales the characteristic directions converge very 
rapidly [7]. This is because the characteristic directions are local quantities, 
and they will typically converge after the fiducial trajectory has explored a 
portion of their neighbourhood. 

In this paper we show that there are constraints on the derivatives of FTLEs 
and of characteristic directions. The derivatives are taken with respect to the 
initial starting coordinate of the fiducial trajectory — the Lagrangian coordi- 
nate — and so are called Lagrangian derivatives. The differential constraints 
are asymptotic, in the sense that they are satisfied with quasi-exponential 
accuracy in time. 

For two- and three-dimensional dynamical systems, these constraints were first 
derived by Tang and Boozer [8] and Thiffeault and Boozer [9]. The method 
used there was to transform the Euclidean metric defined on the space in which 
the motion takes place (Eulerian coordinates) to the space of initial conditions 
(Lagrangian coordinates, or fluid labels). The metric in Lagrangian coordi- 
nates exhibits quasi-exponential blowup, but its associated Riemann curva- 
ture tensor must vanish identically, since the underlying space is flat. Thus, 
by balancing the various terms in the Riemann curvature tensor, it was shown 
that the curvature can only vanish identically if the FTLEs and characteris- 
tic directions obey asymptotic differential constraints. These constraints were 
also confirmed numerically for discrete maps [8,10] and for flows [9,11]. 

Here we propose a different approach, with an emphasis on flows. By per- 
forming a singular value decomposition (SVD) of the tangent map of the 
flow (the Jacobian of the transformation from Eulerian to Lagrangian coordi- 
nates), the quasi-exponentially growing terms can be evolved independently 
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from each other, so that their different growth rates do not pose a problem 
numerically. This amounts to a reorthonormalisation scheme, where the domi- 
nant eigenvectors are subtracted from the subdominant ones [12,7,13]. Greene 
and Kim [14,15] devised a numerical scheme to evolve the components of the 
SVD continuously. Thiffeault [11] extended their method to compute the La- 
grangian derivatives of FTLEs and characteristic eigenvectors. This amounts 
to a numerical scheme for calculating the Hessian (the tensor of second partial 
derivatives) of the transformation from Eulerian (fixed or laboratory coordi- 
nates) to Lagrangian (fluid element labels). Related work on maps was done 
by Dressier and Framer [16] and Taylor [17]. Whereas the Jacobian contains 
information on the deformation of an infinitesimal ball of initial conditions 
into an ellipsoid, the Hessian describes higher-order (nonlinear) deformations. 

The constraints are obtained in two steps. First, the time- asymptotic be- 
haviour of Lagrangian derivatives is derived by analysing the equations for the 
components of the SVD, in a manner analogous to Ref. [11] but allowing for a 
nontrivial (i.e., non- Euclidean) metric. Second, the symmetry of the Hessian 
is used, and Eulerian quantities are replaced by their asymptotic form. What 
is left are expressions generalising the constraints of Ref. [9] to an arbitrary 
number of dimensions and choice of metric. The new method of derivation also 
predicts the convergence rate of the constraints, something that was missing 
from earlier analyses. 

Note that even though it is well-known that the Lyapunov exponents are inde- 
pendent of the particular positive-definite metric chosen to measure them [2,15], 
this is not true of the finite-time exponents and the characteristic directions. 
The details of the metric may affect the short- to moderate-time behaviour of 
the FTLEs greatly, and the characteristic directions are different even in the 
infinite-time limit. Hence, the derivation of the present results for the case of 
an arbitrary metric is relevant from a mathematical standpoint. 

From a physical perspective, metrics that differ from the ordinary Euclidean 
case often arise. Examples are the advection-diffusion equation with an an- 
isotropic, inhomogeneous diffusion tensor, and two-dimensional flow on a de- 
formable membrane [18,19]. Even when the underlying space is Euclidean, 
there are many cases where it is advantageous to use coordinates better suited 
to the problem at hand, which lead to a nontrivial metric (spherical, toroidal, 
rotating, etc.). Hence, there is also physical motivation in keeping the deriva- 
tion of the constraints as general as possible by considering an arbitrary, time- 
and space-dependent metric. 

The outline of this paper is as follows. In Section 2 we review the mathemat- 
ical description of a continuous dynamical system (flow) on a manifold with 
curvature. We define the Jacobian and the Hessian of the flow, and describe 
the metric on the Eulerian and Lagrangian tangent spaces of the manifold. 
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We also choose a convenient covariant time derivative. In Section 3 we define 
the orthonormal frames that will allow the separation of the different scales 
of the system along characteristic directions. This leads to the singular value 
decomposition and the time evolution equations for the orthonormal frames. 
From these we can define the finite-time Lyapunov exponents and characteris- 
tic directions. We devote Section 4 to showing that the nontrivial metric does 
not significantly affect the asymptotic forms of the Lagrangian derivatives de- 
rived in Ref. [11], allowing us to make use of the results therein. Section 5 
contains the main results of the paper: we derive two types of constraints, 
and give the number of constraints of each type as a function of dimension 
and number of contracting directions (negative Lyapunov exponents). In Sec- 
tion 6, we make a connexion to the previous method (the curvature method of 
Refs. [8,9]) of finding constraints. Finally, we offer a summary and discussion 
of the implications of the results of the paper in Section 7. 



2 Chaotic Dynamics in a Curved Space 

In this section we introduce notation and review a number of concepts nec- 
essary for the remainder of the paper. In Section 2.1 we discuss a general 
dynamical system defined on a manifold. In Section 2.2 we equip that mani- 
fold with a metric tensor. To assist in writing equations in manifestly covariant 
form, in Section 2.3 we introduce a covariant form of the time derivative. 

2. 1 Dynamical Systems on a Manifold 

We consider the dynamical system on a smooth compact manifold XL 



where the overdot indicates a derivative with respect to time. Here x = x(t) e 
XL is a smooth function of t G M, and v : It — > TU is a smooth vector field, 
with T XL the tangent bundle of U. We say that v generates a flow </?* : XL — > XL, 
where (p f (a) is a diffeomorphism defined for all a in XL and tel. The flow tp* 
exists, and it satisfies Eq. (2.1) if we set x = tp t (a) [20, p. 304]. We choose ip f 
such that <p°(a) = a, defining a as the initial condition of a trajectory. By 
analogy with fluids, the x are called Eulerian coordinates and the a Lagrangian 
coordinates. 

The derivative (also known as the tangent map, Jacobian, or push-forward) 
at a of </?* is denoted p\ a : T a tl — > T^t^XL, and satisfies the equation 





<fila = 



Dv(t,<p\a))-<p\ 




= Id, 



(2.2) 
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obtained by differentiating Eq. (2.1). We use a dot, as on the right-hand side 
of Eq. (2.2), to indicate the action of a linear map. In a particular chart, 
Eq. (2.2) can be written in coordinate form as 

= ^foOV = (2-3) 

where dv l /dx k is evaluated at x — </?*(a), and we are using the Einstein con- 
vention of summing over repeated indices from 1 to n, where n is the dimension 
of U. 

Equation (2.3) introduces a convention we shall use for the remainder of the 
paper: the indices i,j,k,£ denote vectors in Eulerian space (T X U), while the 
indices p, q, r denote vectors in Lagrangian space (T a U) . Since y?* a is a map 
from T a U to T^U, we write its components as ((pi a ) k q - 

If we imagine an infinitesimally small "ball" of points distributed around the 
initial condition a at t — 0, then the tangent map (fl a captures the deformation 
of that ball into an ellipsoid under the action of the flow </?* (see Ref. [13] for a 
more detailed description). For a chaotic flow, at least one axis of this ellipsoid 
grows exponentially in time, the growth rate being characterised by the largest 
Lyapunov exponent [2,1] (for finite times the growth rate is only roughly ex- 
ponential, or quasi-exponential, but from now on we drop the quasi- prefix). 
This renders the direct numerical evaluation of Eq. (2.3) for moderately large 
times all but impossible due to roundoff error: all the columns of ((p{ a Y q be- 
come aligned with the eigendirection of fastest stretching. A suitable matrix 
decomposition is needed to separate the growth of the different axes of the 
ellipsoid such that the fastest expanding direction does not overwhelm the oth- 
ers, that is, a reorthonormalisation scheme [12,7,13]. In Section 3.1 we achieve 
this using the singular value decomposition (SVD), in a manner analogous to 
Greene & Kim [14,15]. 

In this paper, as in earlier work [11], we are interested in deformations of a ball 
of initial conditions beyond ellipsoidal, so we want to obtain higher derivatives 
of </?*. The second derivative of ip l may be regarded as a symmetric bilinear 
form : T a U x T a U -> T^ (a) U, satisfying 

tpl a = Dv(t, ip\a)) ■ vl a + D 2 v(t, ^(a)M a ® y^J , v?L = 0. (2.4) 

(See Dressier and Farmer [16] for a similar derivation for maps.) Writing out 
the components on a particular chart explicitly, Eq. (2.4) has the form 

(v**a) Pq + (vO P (v?*J , , (y?**j; g = o. (2.5) 

We call the quadratic form ipl^ the Hessian of </?*. The growth rate of the Hes- 
sian is characterised by generalised or higher-order Lyapunov exponents [16,17,11]. 
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These are typically larger in magnitude than the Lyapunov exponents asso- 
ciated with the growth of the tangent map </?* a (the Jacobian). This means 
that numerically integrating Eq. (2.5) directly is even more problematic than 
for the tangent map. In Ref. [11] extensions of the SVD and QR (orthogonal- 
triangular) decomposition methods were introduced to resolve this numerical 
difficulty This allowed the precise verification of the differential constraints 
predicted in earlier work [8,9]. 

2.2 The Metric 

To measure the norm of vectors, we now introduce a Riemannian metric (•, -) x : 
T X U x T X U — > 1RL This induces a corresponding metric (•, •)* on T a U defined 
by 

(X,YY a :=( V i a .X, V i a .Y) x , (2.6) 

for all X, Y e T a U. At t — the two metrics coincide. The reason why we 
introduce a metric on T^U and from it obtain the metric on T a 1i, and not the 
other way around, is best understood in terms of the fluid analogy. We wish to 
measure the growth of the vector X as it is "dragged" by the flow. Initially, the 
norm of X is (X, X)]/ 2 , and at a later time t its norm is (</?* a • X , </?* a • X)]J 2 . 
The latter is then used in Eq. (2.6) to define a time-dependent metric directly 
on T a lX, with no explicit reference to </?*. Thus, it is the metric on T^U that 
is given. 

In component form, we have (X , Y) l a = g pq (t, a) X p Y q , with 

9 Pq (t, a) = (tpi a y p h tJ (t, ip\a)) (cpl a y q , (2.7) 

where hij(t, <^*(a)) are the components of the metric on T X U. In Euclidean 
space, we have hij(t, y?*(a)) = Sij, but we shall retain the more general form. 
Both g pq and hij are positive-definite symmetric matrices, since we chose a 
proper (positive-definite) Riemannian metric. 

2.3 A Covariant Time Derivative 

To ensure that we are always dealing with tensors — and thus guarantee the 
covariance of all expressions — it is convenient to introduce a covariant version 
of the time-derivative operator [21]. In the Eulerian basis, the components of 
this operator D acting on a vector X e T X U are 

M'-r + r^iV + ffc^, (2.8) 
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where the Riemann-Christoffel connexions are defined as 



ifc ._ 1 

ij ' 2 



(2.9) 



and 




(2.10) 



dt 



We use the symbol Vj to denote a covariant derivative with respect to x j , 



Without the 7 term, Eq. (2.8) is the usual definition of covariant differentiation 
along a curve [22], in our case the curve being the trajectory of x with tangent 
vector v. The extra term allows the metric to depend explicitly on time whilst 
preserving the compatibility property Dhij = 0. Note that T, 7, and v do not 
transform as tensors, but their combination defined by Eq. (2.8) does [21]. 



3 Two Orthonormal Frames 

In Section 2.1, we mentioned that direct numerical integration of the evolu- 
tion equation (2.3) for the tangent map ipl a is impractical due to the dom- 
inance of the eigenvector associated with the direction of fastest stretching. 
In the present section we use an appropriate matrix decomposition to sep- 
arates the vastly different scales contained in (p\ a (Section 3.1), as done by 
Greene & Kim [14,15], but our treatment includes a metric with explicit time- 
dependence. We aim to write the equations of motion for the matrices of the 
decomposition in manifestly covariant form (Section 3.2), using the covariant 
time derivative of Section 2.3. In addition to clarifying the role of a nontrivial 
metric in chaotic dynamics, the covariant form allows a more straightforward 
derivation of the equations of motion for the Lagrangian derivatives, as we 
shall see in Section 4. Finally, in Section 3.3 the finite-time Lyapunov expo- 
nents and characteristic directions are defined. 

3.1 Decomposition of <pl a into Orthonormal Bases 

The tangent map </?* a is the Jacobian dx/da of the transformation from La- 
grangian coordinates (a) to Eulerian coordinates (x). Hence, it tells us how an 
set of orthonormal vectors (with respect to the metric h) at a is transformed 
to a set of vectors at x. This last set is not orthonormal in general. 




(2.11) 
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We can however decompose tp\ a into an outer product of orthonormal bases, 

W^EW'W.r (3.1) 

<T,T 

where r] crT = 5 aT , since we have a positive-definite metric, but we write it 
as i] following standard notation [22]. Note that for Greek indices we will not 
use the Einstein sum convention; all sums will be written out explicitly. By 
"orthonormal," we mean that the bases w a and e CT satisfy 

^(M)=EK) ! W 3 T, g pq (t,a)=J2(e,) P (er) q V aT - (3.2) 

<7,T <T,T 

The orthonormal basis vectors e a depend on a and t, and the basis vectors w T 
depend on x and t; they diagonalise g and h, respectively. The basis {e^} 
is a noncoordinate basis, i.e., it does not correspond to the natural tangent 
vectors {d/da q } of the coordinate system a, or of any other coordinate system. 
The same is true of the basis {w a }. It follows from Eq. (3.2) that 

(w T Y (w a ) t = rj Ta , (e T ) q (e a ) q = r] Ta . (3.3) 



In the appendix we show that, using the standard singular value decomposition 
(SVD), we can always write 

M\ = Y,U\F^V qT) (3.4) 

<7,T 

where the matrix M denotes the components of ^>\ a in a particular chart. The 
notation for the components of the decomposition is that of Refs. [14,15,11]. 
The matrix F is diagonal, and the matrices U and V satisfy the orthogonality 
relations, 

hi, U\ U\ = r) aT , rT V pu V qT = VaT . (3.5) 

Note that U is orthogonal with respect to the metric on T X U, but V is or- 
thogonal with respect to rf q , which is not the metric on T a U. The reason for 
this is that we want to avoid including in V the exponentially growing terms 
contained in g. 

The decomposition (3.4) is unique up to permutations of rows and columns. 
The diagonal elements F aa —. A a are called the singular values of M. Re- 
quiring that the singular values be ordered decreasing in size makes the de- 
composition unique (for nondegenerate eigenvalues). We refer the reader to 
Refs. [13,11] for a discussion of the geometrical significance of the SVD in 
dynamical systems. 
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From Eqs. (3.1) and (3.4), we can make the identification 



(e ff ) p = V pp F p a , (e a y = V% {F~ 1 ) p a , (3.6) 

where we used r] pu as a "metric" to raise and lower Greek indices, and we 
defined V p p := i] pq V qp to abridge the notation. Note that V p p is a tensor, but 
is not equal to V p p = g pq V qp because r] pq is not the metric on T a lX; hence the 
need for a different symbol. 

Analogously to {e a }, we can take the basis {w a } to be 

(w a y = U\, (w a )i = U iT . (3.7) 

With the definitions (3.6) and (3.7) it is easy to check that {e CT } and {%} 
satisfy the orthonormality conditions (3.2) and (3.3). 



3.2 Time Evolution 



Greene & Kim [14] derived equations of motion for U, F, and V for a met- 
ric without explicit time-dependence. Our derivation parallels theirs, except 
for the use of the covariant time derivative of Section 2.3 to allow for time- 
dependence of the metric. 

We differentiate the decomposition (3.4) using the covariant derivative (2.8), 
DM\ = = J2 l® 17 ^ paT V qt + U\ T>F aT V qT + U\ F <JT T>V qT ] , (3.8) 

<T,T 

and then contract with U ip V q u , 

Y, [Ui» VU\ F aT Vru + r)^ F UT t) tv + Vpa F aT V\ DV qr ] = 0, (3.9) 

<7,T 

where because F UT is a scalar we replaced the covariant time derivative by an 
ordinary one. To simplify the notation, we use the diagonality of F, and raise 
and lower Greek indices with 77, yielding 

U ip VU\ F\ + F pv + F% V q v VV m = 0. (3. 10) 

The scalar U ip T>U\ is antisymmetric in \x and v, but not V q v DV qil , be- 
cause r Dr] pq 7^ since rf q is not the metric. However, the quantity V q v V qp = 
V pq Vpv Vqu is antisymmetric in ji and u, so we expand the covariant derivative 
of V, 

U ip DU\ F\ + F pu + F\ V q u V m - F\ G pv = 0, (3.11) 
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where 



G^:=U\U\[V j v i + v j }, (3.12) 

and we used the transformation law from Lagrangian to Eulerian coordinates 
imposed by the covariance of D [21]. The matrix G is a covariant generalisation 
of the velocity gradient tensor, expressed in the U basis. The 7 term denotes 
straining due to the time-dependence of the metric. 

We can now solve Eq. (3.11) for the various time derivatives by using the 
antisymmetry of T>U e u and V q ^ V qu . We obtain finally 

F% = G%F%, (3.13) 
V* = - ^^j!!^(/ )2 ; for „ + v; (3.14) 

VQ ^ = ' { F^f-{F%f A ^ ior^v; (3.15) 

where := G^ + G Ufl . These equations are identical in form to those of 
Greene & Kim [14,15], except for the definition of G and the use of the opera- 
tor D instead of the time derivative in the U equation. Note that the matrix A 
is the rate-of-strain tensor (up to a possible factor of two, depending on the 
convention used), expressed in the basis U . 

The ordinary time derivative is used for F because it is a scalar. For U, we need 
to use the modified version D of the time derivative to take into account the 
nontrivial, possibly time-dependent metric. The ordinary derivative is used 
for the V equation, even though it is a vector, because we have effectively 
rescaled the metric g pq to give r] pq . The reason for doing so is that for a chaotic 
flow the elements of g pq grow exponentially, rendering the metric difficult to 
use directly. Instead, we have absorbed the exponential growth into the F^, 
which are the coefficients of expansion of the flow [2,1]. 

The techniques used in this section also apply to the QR method for obtain- 
ing the Lyapunov exponents and characteristic directions [13,7]. One simply 
replaces the ordinary time derivative of the Q equation by a covariant deriva- 
tive (analogous to the U), and the time derivatives of the elements of R are 
unchanged (analogous to F and V). The dv/dx term must be modified to 
use the covariant spatial derivative, and a time derivative of the metric must 
be included, exactly as in Eq. (3.12). The QR method has the advantage of 
having no singularity in its equations of motion, as opposed to Eqs. (3.14) 
and (3.15) which are singular for = A u . However, the type of asymptotic 
analysis we do here is better done with the SVD method, as it was in Ref. [11]. 
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3.3 Lyapunov Exponents and Characteristic Directions 



The relation (3.6) between the orthonormal frame {e a } and the SVD matrices 
allows the definition of the finite-time Lyapunov exponents 

\*(t,a) := - t logA^a), (3.16) 

where the A CT := F aa are called the coefficients of expansion [2]. The A CT 
give the instantaneous relative growth of the principal axes of an infinites- 
imal ellipsoid moving with the flow. Taking the limit t — > oo in Eq. (3.16) 
gives the Lyapunov exponents A£°, which are independent of a and t in a 
given ergodic region [2] (for almost all initial conditions); they converge very 
slowly [7]. Associated with the coefficients of expansion are the characteristic 
directions (e a ) q = V qa , which converge exponentially rapidly to their time- 
asymptotic value e^°(a) [7,11] (for nondegenerate Lyapunov exponents). They 
give the directions of stretching of the axes of the ellipsoid in Lagrangian 
space. 

Following Goldhirsch et al. [7], we shall assume that the Lyapunov expo- 
nents are nondegenerate and ordered such that A (7 _ 1 > A a . After allowing 
some time for chaotic behaviour to set in, we have that A CT ^> A K for a < k. Use 
will be made of this ordering in the next section to estimate the asymptotic 
behaviour of Lagrangian derivatives. 



4 Lagrangian Derivatives 



In Thiffeault [11], evolution equations for the Lagrangian derivatives of e (V), 
w (U), and A (F) were derived. In this section we show how these equations 
must be modified to account for a nontrivial metric. We then demonstrate that 
the asymptotic form for the Lagrangian derivatives, also obtained in Ref. [11], 
applies to the curved case. 

We define the quantities 

* w := V p K ^-logF\, := V\U\V p U iv , (4.1) 



e K „ := V\ V\ ^ (4.2) 

which are the Lagrangian derivatives of A u , Ui U , and V qu expressed in the 
orthonormal frames {w a } and {e a }. Integration by parts leads to the sym- 
metries Q Kflu = —Q KUl i and Kflu = —Q KV fx- We are using an ordinary spatial 
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derivative — as opposed to a covariant derivative — in the definition of 9 be- 
cause we have effectively rescaled the metric g pq to give r) pq (Section 2). Using 
the covariant derivative would destroy the antisymmetry of in its last two 
indices. 

Equations of motion for $, and G can be derived in a manner analogous to 
Ref. [11], by taking Lagrangian derivatives of the equations of motion (3.13)- 
(3.15). However, care must be taken when commuting covariant derivatives: 
since the metric h is nontrivial, curvature terms can arise. Moreover, a form 
of time-curvature associated with the nontrivial time-dependence of h is also 
present [21]. The evolution equation for $ is 



= E (v T vy K <s>^ + E (u T m)\ - E (u T my u $ 



Kfia 



+ V\ -^- q (U T VU), u - K QT R^ Ka v° + S, UK ) (4.3) 

where the matrices (U T 1)U) and (V T V) denote the right-hand side of Eqs. 
(3.14) and (3.15), respectively. The last term in Eq. (4.3) appears only when h 
has nontrivial spatial and temporal dependence. It contains the curvature R, 
defined by [22] 

E v ° = U\ U\ U\ [V k V, Vj - V, V k Vj ] , (4.4) 



and the time-curvature S [21], defined by 

S^ K := U\ U\ W K [V fc 7ii - V i7jfc ] • (4.5) 
Both R^ UKa and S^ VK are antisymmetric in ji and v. 
Equations of motion can also be found for \l> and 6, 

= E (v T v)\ * w + E A ^ + A - x — > ( 4 - 6 ) 



= E (^vT, + E (^vT„ - E (^vT„ 

cr a cr 



+ V^-^(V T VU, (4.7) 

where 

:= U\ U\ U e , V, [V, v k + 7M ] . (4.8) 

Equations (4.6) and (4.7) are identical in form to their Euclidean versions of 
Ref. [21], with the proviso that the X in that paper must be replaced by X 
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of Eq. (4.8). In its new form, which allows for a nontrivial metric, X UKfl is no 
longer symmetric in k and rather, it satisfies the commutation relation 

■^VKfl -^VflK ^ ' R^IVKIT V . (4-9) 



Without redoing the entire derivation, we have thus shown that the evolution 
equations derived in Ref. [11] can be easily modified to allow for a nontrivial 
metric h. Moreover, it is straightforward to show that these terms do not 
modify the asymptotic behaviour of the Lagrangian derivatives derived in 
Ref. [11]. The reason is that if we assume that the motion takes place on some 
bounded region of phase space, then the curvatures R and S must be bounded. 

We now quote the result of Ref. [11], which also applies to our system: for 
t ^> 1, by which we mean that the dynamical system has evolved long enough 
for the quantities A M to have reached a regime of quasi-exponential behaviour, 
we have that the Lagrangian derivatives defined by Eqs. (4.1) and (4.2) evolve 
asymptotically as 

= max (A K , 7/1 „) $ Kflu , (4.10) 
= max (A K , 7kk , v) § KV + , (4.11) 
9 K/li , = max (7 MJ ,A K , 7 KK , 7 M/1 , K ^ + 6^ , (4.12) 



where 



7^ : = \ 



A,/A u , ix>v- (4.13) 



The quantities with tildes have algebraic (nonexponential) behaviour, and the 
quantities with 00 superscripts are constants. The 7's are defined such that 
asymptotically they are decreasing exponentially, for nondegenerate A. See 
Ref. [11] for a more detailed discussion of the character of this asymptotic 
behaviour. 



5 Constraints 



The Hessian y?'^, introduced in Section 2, is the tensor of second derivatives 
of y?*(a); it characterises deformations of fluid elements beyond ellipsoidal. 
Because partial derivatives commute, the Hessian is symmetric in the corre- 
sponding indices. In Ref. [11] we showed that this symmetry of the Hessian 
implies that the Lagrangian derivatives of Section 4 are not all independent. 
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However, the relationship between the derivatives is very singular for chaotic 
flows. We show in this section that the relations can be used to obtain con- 
straints on the derivatives, recovering and extending previous results [8,10,9] 
deduced by different means. 



5. 1 Symmetry of the Hessian 

We define the projection of the Hessian (y?** a )p g = d 2 x e /daPdcfl onto the U 
and V bases as 

K K , U := U tK (<pl a Y pq V% V\ = K KV , . (5.1) 
Since (<p{ a Y q = dx e /da g , we can write 

daPdai daP ' 1 ' 1 

and using the decomposition (3.4) for M % q = (fl a Y q , we find 

K KfMU = A K r] UK + A K Q^ VK + A u . (5.3) 

The Hessian is symmetric in its lower indices, so we can equally well write 

K K[XV = A K ^ VK r]^ + A K 6^ K + A M $ VKIi , (5.4) 

where we interchanged /i and v in (5.3). Equating (5.3) and (5.4), we find the 
relations 

A M + V^v^) = K <$>^ v , /i^, (5.5) 

A K (Q^uk - Qu^k) = A u $n UK - A M , /I, v, k differ. (5.6) 

Equations (5.5) and (5.6) can be used to solve for the G's in terms of the $'s 
and \l/'s. However, for chaotic flows the exponential behaviour of the A's ren- 
ders this solution highly singular [11]. Thus, in practise it is better to evolve 
all the Lagrangian derivatives. 

Rather than solving for the 0, if the flow is chaotic the relations (5.5) and (5.6) 
can be put to good use in another manner. The Lagrangian derivatives of U, 
as contained in <£>, are not quantities of great interest to us. They contain in- 
formation about the sensitive dependence on initial conditions of the absolute 
orientation of fluid elements in Eulerian space. This information is not nec- 
essary for solving problems in Lagrangian coordinates, and is too sensitive to 
initial conditions to be of use anyhow. We thus substitute the time-asymptotic 
form of <3>, given by Eq. (4.10), in the right-hand side of Eqs. (5.5) and (5.6), 



14 



yielding 

Q^u + Vw®^ = max , ^ 7^ , (5.7) 
Q^vk - Qvhk = ~t~ max (A M , n/ VK ) $ M „ K - ^ max (A„ , 7 MK ) , (5.8) 



where, as before, /i, v, and k differ. 
5.2 Type I Constraints 

For /j, < v, Eq. (5.7) can be written 

+ = max (a„ , 7^) , n < v, (5.9) 

where we have set rj^ = 1 for a proper Riemannian metric. If the index z/ cor- 
responds to a contracting direction (Aj, <C 1), the right-hand side of Eq. (5.9) 
goes to zero exponentially fast, at a rate A u or 7^,, whichever is slowest. In 
that case Eq. (5.9) is a constraint implying that for large t we have 

(6^ + — ► 0, A, « A M , A, < 1. (5.10) 

We refer to (5.10) as Type I constraints. The total number of such Type I 
constraints is 

Ai = n s [n- 5K + I)] , (5.11) 

where n is the dimension of the space and n s is the number of contracting 
directions (i.e., the number of negative Lyapunov exponents) possessed by 
the flow in a particular chaotic region. Table 1 gives the number of Type I 
constraints, Ni, as a function of n and n s . 

In two dimensions, we typically have one contracting direction, so there is 
a single Type I constraint. This is the same constraint that was derived in 
Refs. [8,9]. 

In three dimensions, for an autonomous flow, we typically also have one con- 
tracting direction. There are then two Type I constraint. These constraints 
correspond to those derived in Ref. [9]. 

An interesting special case of the Type I constraints is obtained by set- 
ting v = n in Eq. (5.9), and then summing over \i < n. We obtain 

V, {e n y - (e n y V q log A n ~ max (A„ , 7 * n ) - 0. (5.12) 
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Table 1 

The total number of Type I constraints for low-dimensional systems, as given by 
Eq. (5.11). The rows denote n, the columns n s < n. 



n 


1 


2 


3 


4 


5 


6 


7 


1 

















2 


1 


1 












3 


2 


3 


3 










4 


3 


5 


6 


6 








5 


4 


7 


9 


10 


10 






6 


5 


9 


12 


14 


15 


15 




7 


6 


11 


15 


18 


20 


21 


21 



This constraint was discovered numerically and used by Tang and Boozer [8,10,23,24] 
and derived in three dimensions by Thiffeault and Boozer [9]. The present 
method not only gives the constraint in a very direct manner for any di- 
mension n, but it also provides us with its asymptotic convergence rate, as 
determined by the right-hand side of (5.12). 



5.3 Type II Constraints 



Equation (5.14) implies that 

- ~ max — , — 7^ , — 7^ j . (5. 13) 

We are interested in finding constraints analogous to the Type I constraints 
of Section 5.2. It is clear that unless both fi and v are greater than k, the 
right-hand side of Eq. (5.13) is of order unity or greater, and so does not go 
to zero. We can assume without loss of generality that /i < v, so that 

Q^uk ~ ~ lp,K max (A„ , 7 MK ) , v > fi > «, (5.14) 

where we have used r ) VR <^ 7 MK . Whether or not Eq. (5.14) is a constraint 
depends on the specific behaviour of A M K v j K K . Clearly, we have a constraint 
if K v 1, since 7 MK 1. This provides a lower bound on the number Nn of 
Type II constraints; by choosing v from the n s contracting directions, and 
summing over the remaining k > /i > v, we obtain 

N u > \n s [n 2 - (n s + 2) (n - \{n s + 1))] . (5.15) 
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Table 2 

Lower bound on the number of Type II constraints, as given by Eq. (5.15). The 
rows denote n, the columns n s < n. 



n 


1 


2 


3 


4 


5 


6 


7 


1 

















2 


















3 


1 


1 


1 










4 


3 


4 


4 


4 








5 


6 


9 


10 


10 


10 






6 


10 


16 


19 


20 


20 


20 




7 


15 


25 


31 


34 


35 


35 


35 



But even if A u ^> 1 we can have a constraint, as long as A u / ~f^ lK <C 1. This 
depends on the particular problem at hand; hence, Eq. (5.15) is only a lower 
bound, but a fairly tight one for low dimensions. Table 2 enumerates the 
minimum number of Type II constraints as a function of n and n s . 

Note that when k, //, and v differ we can write 

®fu/K ~~ ®i//j,k = ~ (^k)<j ? ^u] q , (5.16) 
where the Lie bracket is 

& , ^ := {e,f A _ (e ^)P ^_ (6/t)9 . ( 5 .17) 

The Type II constraints are thus forcing certain Lie brackets of the character- 
istic directions e a to vanish asymptotically. The geometrical implications of 
this, and perhaps a connexion to the Frobenius theorem and the existence of 
submanifolds, remains to be explored. 



6 Curvature 



We now make contact with the earlier results of Refs. [8,9], where the con- 
straints were derived for two and three dimensional flows by examining the 
form of the Riemann curvature tensor associated with the metric g. We list the 
advantages of the present method. The only known disadvantage is that the 
formalism presented here does not apply to maps, only to continuous flows. 
This is because of the lack of an SVD method for maps [13]. The use of the 



17 



QR method, which is defined for maps [1,13], may circumvent that problem, 
but the analysis is more difficult and has not been carried out yet. 



Nontrivial metrics can have curvature, as reflected by the presence of the 
Riemann curvature tensor in Section 4. A straightforward method of obtaining 
that tensor is through the use of Ricci rotation coefficients [22], 

u>^v := (w K y (w p ) j V, {w v )j = {e K f {e p ) q V p (e„), . (6.1) 

These satisfy the antisymmetry property uj Kpu = —ou Kup , and can be rewritten 
in terms of the $ of Eq. (4.1) as 



CO 



'.-/<" -A" 1 ^. (6.2) 



In terms of the rotation coefficients, the Riemann curvature tensor is [22, p. 51] 

— f] TP l^Kpp ^(jtv — W apll LO KTV + 0J Kpa U) Tpv — L0 apK LO Tpv \ . (6.3) 

If we use the relations (5.5) and (5.6) to solve for $ in terms of 6 and \P, we 
can rewrite Eq. (6.2) as 

^Kpu = 2t\ \ \ {^-^(©i/K/i — Qnup) + Aj,(0 K/1 i, — Q^ KV ) — A K (&^ UK — @^ K )| 

+ 1 -r]^up- 1 -r]u^pu (6.4) 

The form of the curvature obtained by inserting Eq. (6.4) into (6.3) is es- 
sentially the one obtained in three dimensions in Ref. [9]. The curvature was 
calculated directly from the Christoffel symbols in Ref. [8]. The constraints 
were then deduced by imposing the finiteness of the curvature tensor: some 
terms in Eq. (6.4) would appear to grow exponentially, so their coefficient 
must go to zero to maintain the curvature finite. In this manner, the Type I 
and Type II constraints were derived in two and three dimensions, backed by 
numerical evidence [8,9,11]. 

The approach used in this paper to derive the constraints is advantageous in 
several ways: (i) It is valid in any number of dimensions; (ii) It avoids using the 
curvature; (iii) There are no assumptions about the growth rate of individual 
terms in the curvature [9]; (iv) The convergence rate of the constraints is 
rigorously obtained [Eqs. (5.9) and (5.14)]; (v) The number of constraints can 
be predicted [Eqs. (5.11) and (5.15)]. The crux of the difference between the 
two approaches is that here we use the variable $ to estimate the asymptotic 
behaviour of the constraints directly, rather than relying on indirect evidence 
from the curvature. Thus, the derivation of the time-asymptotic form of $ is 
essential (Section 4 and Ref. [11]). 
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7 Discussion 



In this paper we have demonstrated that in chaotic flows the Lagrangian 
derivatives of finite-time Lyapunov exponents, and their associated charac- 
teristic directions, satisfy time-asymptotic differential constraints. The con- 
straints are valid for any metric tensor, and are satisfied with quasi-exponential 
(that is, exponential in the long-time limit) accuracy in time. The method 
employed to derive the constraints made use of earlier work on the asymp- 
totic form of Lagrangian derivatives [11]. The new results generalise previous 
work [8,9] by providing constraints in any number of dimensions, together with 
their convergence rate. We now give a few examples of possible applications 
of the constraints, and directions for future research. 

In Ref. [8] it was demonstrated that for the isotropic, homogeneous advection- 
diffusion equation the diffusion proceeds almost exclusively along the charac- 
teristic direction associated with the smallest Lyapunov exponents. This was 
shown by transforming the equation to Lagrangian coordinates, thereby elimi- 
nating the advection term. The resulting equation was an anisotropic diffusion 
equation, where the diffusion tensor was proportional to g -1 , the inverse of the 
metric tensor in Lagrangian space. Because the inverse of the metric appears, 
the dominant eigendirection is that associated with the smallest Lyapunov 
exponent. Physically, this means that diffusion occurs in the direction along 
which large gradients are created by the stretching effect of the flow. 

If the diffusion tensor in the advection-diffusion equation is anisotropic and 
inhomogeneous, then the direction of fastest diffusion is associated with the 
characteristic direction calculated using the diffusion tensor itself as a metric. 
Furthermore, in the large Peclet number limit (small diffusivity) , the Type 
I constraints derived in Section 5.2 can actually be used to transform the 
advection-diffusion equation into a scalar one-dimensional diffusion equation. 
This use of the constraints will be addressed in subsequent work [19]. Other 
applications of the Type I constraint can be found in Tang and Boozer [23,24]: 
they demonstrate that regions of high curvature of the integral curves of e n 
lead to locally small finite-time Lyapunov exponents, hindering the chaotic 
enhancement of diffusion. 

In three dimensions with an Euclidean metric, the Type II constraint (5.14) 
can be written 

ei • V x ei ~ 712 max (A 3 , 712) 0, (7.1) 

where the curl is taken with respect to Lagrangian coordinates, a. Equa- 
tion (7.1) implies that, asymptotically, e.\ is of the form 

ex = W>/ IIVVU , (7.2) 
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where ||-|| is the Euclidean norm, and ^ is a pseudopotential function. Equa- 
tion (7.2) is valid locally. Hence, the field of unstable directions (associated 
with the largest Lyapunov exponent) has a simple asymptotic structure. The 
implications of that structure remain to be investigated. An indication of how 
it might prove useful was investigated in Ref. [9], where it was shown that 
the Type II constraints could have an impact on energy dissipation in the fast 
kinematic dynamo. 
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A Singular Value Decomposition with a Metric 



The decomposition given by Eq. (3.4) is denned for a space with a nontriv- 
ial metric, but the numerical algorithms for obtaining the SVD are based 
on "true" orthogonal matrices (that is, matrices orthogonal with respect to 
the Euclidean metric). In this appendix we show that we can always find 
the bases {e CT } and {w a } using the ordinary singular value decomposition of 
matrices with a Euclidean metric. 

Since the matrix of components of the metric (•, -) x on T X U is symmetric, 
it can be diagonalised by an orthogonal transformation. This transformation 
can then be composed with a rescaling to make the diagonal elements equal 
to unity (a proper Riemannian metric has positive eigenvalues), yielding 

hij = k iTji'ji k? j, (A-l) 

where k 1 ' \(t, x) is the coordinate transformation, and r^y = S^f for a positive- 
definite metric. We are using primed indices to denote this new basis on T X U 
where h is diagonal. The metric g of Eq. (2.7) can then be written 

gp^^taYp^iVi'fk^UJ,. (A.2) 

Define the matrix M with components 

M% := k\ (ri a )\ , (A.3) 

so that Eq. (A.2) can be written 

g pq = M i ' p r H , f M j ' q . (A.4) 
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Since r^y = (5^, the metric is now in a form g = M T M, suggesting the use of 
the usual SVD technique to eliminate the Eulerian information in g. We have 
absorbed the non-Euclidean metric into M, so that even though the space is 
not Euclidean the decomposition carries through in the same manner as in 
Refs. [14,11]. 

Proceeding with the singular value decomposition, we write 

M' q = Q l ' a F° r V qr , (A.5) 

where Q and V are ordinary orthogonal matrices (i.e., with respect to the 
Euclidean metric) and F is diagonal. 

The basis {e CT } is given by V and F as in Eq. (3.6). The basis {w a } is given 
in terms of Q and k by 

(w a y = U\ = (AT 1 )*, Q l ' CT , K), = U ia = m,? V'i Q j ' a . (A.6) 

Since the SVD always exists, we have demonstrated that the construction of 
the basis vectors e a and w a satisfying Eqs. (3.1) and (3.2) is possible. 
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